The X-Philes (2nd Revision)/The X-Philes Number 1 (1995).iso/xphiles/hp48hor1/sunclk.src

Text File | 1991-10-19 | 20KB | 739 lines

%%HP: T(3)A(R)F(.);
@ SUNCLK, by James Elliott
@ Directory: SunClk
@ Checksum: # 40750d  Bytes: 8805.5
DIR
  Now
@ Put the current date and time into TM and DT, convert to Julian
@ date/time, and call SunClk, the driver program.
    \<< DATE TIME
'TM' STO 'DT' STO
GTime 1 'STAT' STO
SunClk
    \>>
  Then
@ Create a temporary menu that allows specification of an arbitrary
@ date and time, and display the current settings. Trap errors in case
@ an invalid date or time is entered, so garbage isn't left on the
@ stack.
    \<< { { ">DATE"
@ Menu key for setting the date to view.
      \<< \-> i
        \<< DT \-> d
          \<< i
            IF DUP
DUP IP ==
@ If only the month was entered, default the rest from the previous
@ setting.
            THEN d
FP +
            ELSE
              IF
DUP FP 100 * DUP IP
==
@ Similarly, if the year is missing, provide the default.
              THEN
d 100 * FP 100 / +
              END
            END
'DT' STO
@ Update the date to reflect their changes, fixing things if there is
@ an error.
            IFERR
ShwTm
            THEN
ERRM 3 DISP 2 WAIT
DROP DROP d 'DT'
STO ShwTm
            END
          \>>
        \>>
      \>> } { ">TIME"
@ Menu key for setting the time to view. Very similar to above.
      \<< \-> i
        \<< TM \-> t
          \<< i 'TM'
STO
            IFERR
ShwTm
            THEN
ERRM 3 DISP 2 WAIT
DROP DROP t 'TM'
STO ShwTm
            END
          \>>
        \>>
      \>> } { "A/PM"
@ Toggle AM vs. PM
      \<< TM 12 - DUP
        IF 0 <
        THEN 24 +
        END 'TM'
STO ShwTm
      \>> } { } { } {
"Go"
@ Use the specified time; convert it to Julian time and call SunClk,
@ the driver program, restoring the normal menu.
      \<< 0 MENU
GTime 0 'STAT' STO
SunClk
      \>> } } TMENU
@ Okay, the temporary menu is built; now show what the default time
@ and date are (these are the ones used for the last display.)
ShwTm
    \>>


@ Provide brief instructions. Blank the menu while waiting for
@ keypresses between screens.
  Help
    \<< CLLCD { }
TMENU
"Press NOW for current"
1 DISP
"daylight map, or THEN"
2 DISP
"to pick a day & time."
3 DISP
"Set TZ to your time"
5 DISP
"zone (eg. CDT = -5)."
6 DISP -1 WAIT DROP
CLLCD
"To interrupt, press"
1 DISP
"any key other than"
2 DISP
"ON/ATTN, and it will"
3 DISP
"stop in a few seconds"
4 DISP
"and clean up." 5
DISP 0 WAIT DROP
CLLCD
"'NOW' updates the map"
1 DISP
"continually until you"
2 DISP
"interrupt it, 'THEN'"
3 DISP
"draws a map and ends."
4 DISP 0 MENU 3
FREEZE
    \>>


@ Provide brief credits. Blank the menu while waiting for
@ keypresses between screens.
  About
    \<< CLLCD { }
TMENU
"HP 48SX implementation"
1 DISP
"by James J. Elliott"
2 DISP
"<[e-mail]>"
3 DISP
"Public domain; share"
5 DISP "and enjoy!"
6 DISP -1 WAIT DROP
CLLCD
"Based on a SunView"
1 DISP
"program by:" 2
DISP "John Walker"
4 DISP
"<[e-mail]>"
5 DISP 0 WAIT DROP
CLLCD
"Algorithms found in"
1 DISP
"'Astronomical Formulae"
2 DISP
"for Calculators' by"
3 DISP
"Jean Meeus, Third"
4 DISP
"Edition, 1985." 5
DISP 3 FREEZE 0
MENU
    \>>


@ This variable contains your current time zone. My default is Central
@ Daylight savings Time, or 5 hours earlier than UTC.
  TZ -5

@ This variable determines how detailed (and therefore slow) the
@ computation of the daylight bands' width is. In projecting the
@ sunlight region over the daylight hemisphere, this many different
@ bands are computed, and linear interpolation is used between them up
@ to the resolution of the graphics display. Note that it does not
@ make sense to set this higher than 62, since only 62 values are
@ stored! The lower the number, the chunkier the display, but the
@ faster arbitrary dates can be displayed. It doesn't really affect
@ the speed of ongoing updates, since this computation is only done
@ when the sun's declination drifts beyond a threshold from its
@ current value.
  Res 40

@ This is the threshold referred to above.
  Thres .5

@ This is just a band of black used in drawing the sunlight (with
@ GXOR). It's much faster than using LINE.
  STRIP
GROB 131 1 FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF0

@ Here's the world map.
  World
GROB 126 64 FFFFFFF30000008FFFFFFEFFFFFFFFF3FFFFFF1000CFF18F10F00CF30FFFFFF3FFFF7000008FF7EF18FF38F30CF1EFF3FFFF70000E3FF1EFFFFF002840018FF3E30E1000007CF07EF0CF010FF5000892C8F00000004C78FF340000CF9FFF910208FF744000CC038F1008F30FFFFFFF72389FF0C000E0E3C7800EFFFFFFFFF132F000E18971E1FFF0000FFFFFFFF70083F3C3C72324CFFFF8004FFFFF9FF1F0F370FF8F7E0D1FFF700CF3FFCF1EF1E8F30EFF0F74270FFF70FFF9FFDF3FF7ECF3FFFF3E17020EFFF9FF7D7EFFFFF36EF3FFFF7F73014EFFFB0F03270FFFF10FF3FFFF7F7F00EFFF7006122FCFFF74CFF3FFFF7FFF88FFFE730024CFFFF703FFF3FFFF7EFFFEFFF8700184EFFFF701FFF3FFFFF8FF3EFFFF7883BFFFFFF720FFF3FFFFF3F78FFFFF4F3083FFFFF76EFFF3FFFFF3E09FFFFF1FFF93CFFFFFEFFFF3FFFFF3CE3FFFFF9FFF170CFFF7EFFFF3F7FFFF8E0FFFFFCFFF0F81F870EFFFF3F7EFFFB000FFFFCFFF3E93727EEFFFF3FFFFFF3040EFFFDFFF7EC73766EFFFF3FFFFFFF1CFFFFFCF3F70E7A346CFFFF3FFFFFFEFC0CFFF9F3FF1FF8B06CFFBD3FFFFFFEF008FFFBBFFF7FF8B078FFBD3FFFFFFFF3F3EFF30CFB3FF93019FFFF3FFFFFFFF7FFCFFFFCF89FFF7001FFFF3FFFFFFF33FF0FFFFCE8CFFF7820AFF73FFFFFFFB3FF4CFFFDE4EFFFF00002F73FFFFFFFFBFFF9FFF9F6FFFFF10E90EF3FFFBFFFF3FFFBF9FB76FFFFF708388F3FFFBFFFF7EFF9FFFBF64FFFFF7448FF3FFFFFFFFFAFFDFFFBFC0FFFFFF15EFE3FFFDFFFFF8FFDFFC3300FFFFFFC1CFC3FFFDFFFFFBFFCFFF3BB1DFFFF3EF97C3FFFFFFFFFBF3EFFF3FB9DFFFF9FF3FE3FFFFFFFFFBF9FFFF7F99FFFFF9FF7FF3FFFFFFFFF9F9FFFF7ECFFFFFF1E17FF3FFFFFFFFFDFCFFFFF6EFFFFFF3107FF3FFFFFFFFFD3EFFFFF0FFFFFFF3C93F33FFFFFFFFFD3FFFFFFFFFFFFFFFF18F32FFFFFFFFFD8FFFFFFFFFFFFFFFFFCF32FFFFFFFFFCEFFFFFFFFFFFFFFFFFCF03FFFFFFFFF4EFFFFFFFFFFFFFFFFFFF83FFFFFFFFF4FFFFFFFFFFFFFFFFFFFFE3FFFFFFFFF09FFFFFFFFFFFFFFFFFFFF3FFFFFFFFF0EF9FFFFFFFFFFFFFFFFFF3FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF3FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF3FFFFFFFFFF9FFFFFFFFFFFFFFFFFFFF3FFFFFFFFF7CFFFFFFFFF0E300000EFF3FFFFFFFF70CFFFFF7F00E18FFFFF08F3FFFFFF7040DFFF7000EFFFFFFFFFF3C3FFF1000C08CFF30FFFFFFFFFFFFFF3C3F70C9FFFF70E38FFFFFFFFFFFFFFF3E3F90FFFFFFFF1CFFFFFFFFFFFFFFFF1E3F7EFFFFFFFFFFFFFFFFFFFFFFFFFF1028F9FFFFFFFFFFFFFFFFFFFFFFFFFF7C3FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF3FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF3FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF3FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF3

@ (End of amendments)

@ This reverses the pixels of a region of the line being built up,
@ starting at the pixel in level 2, and ending at the pixel in level
@ 1. It's called by XSpan, below.
  XLine
    \<< \-> s f
      \<< STRIP f 1 +
s - R\->B # 0d 2
\->LIST STRIP GXOR s
2 + R\->B # 0d 2
\->LIST SWAP GXOR
      \>>
    \>>

@ Reverse a region of the line, centered on the pixel in level 1, and
@ twice the length specified in level 2. Handles wrapping off the left
@ edge if needed.
  XSpan
    \<< \-> l n
      \<< l 126 MOD
'l' STO
        IF l n +
126 >
        THEN l 125
XLine 0 l n + 127 -
XLine
        ELSE l l n
+ 1 - XLine
        END
      \>>
    \>>


@ Draw the time at the bottom of the map.
  XTime
    \<< \-> s
      \<< s 1 \->GROB
DUP SIZE DROP \-> o w
        \<< PICT 131
w - 2 / # 59d 2
\->LIST o GXOR
        \>>
      \>>
    \>>


@ Update the sunlit region display if necessary, and the time display
@ if appropriate. Takes the pixel location of noon in level 1.
  Draw
    \<< \-> n
      \<< 0 255 \-> w o
        \<<
@ Is this an update, or the first display? If an update, set flag 6,
@ so that the old lines will be taken into account.
          IF 'OTAB'
VTYPE 0 \>=
          THEN 6 SF
@ Has the sun moved, as far as the resolution of the map is concerned?
@ If so, set flag 7.
            IF n
ONOON \=/
            THEN 7
SF
            ELSE 7
CF
            END
          ELSE 6 CF
7 CF
          END 1 64
@ Loop over the whole map, one horizontal band at a time.
          FOR i
            IF 6
FS?
@ If it's an update, determine what the width was last time.
            THEN
OTAB i DUP SUB NUM
'o' STO
            END
@ Figure out the width of sunlight at this latitude, by looking it up
@ in the width table built by PrjIll.
WTAB i DUP SUB NUM
'w' STO
@ If the sun has moved or the width has changed, we need to update
@ this band.
            IF o w
\=/ 7 FS? OR
            THEN
@ Create a blank line in which the regions which need to change will
@ be marked.
# 131d # 1d BLANK
@ If there was an old sunlit band, we'll erase it.
              IF o
255 <
              THEN
ONOON o - 126 + 126
MOD o 2 * XSpan
              END
@ If there is a new sunlit band, draw it.
              IF w
255 <
              THEN
n w - 126 + 126 MOD
w 2 * XSpan
              END
@ Take the net result of the above two operations and apply it to the
@ actual map in one fast GXOR, so there's no flicker.
PICT SWAP # 0d i 1
- R\->B 2 \->LIST SWAP
GXOR
            END
@ If the user wants to stop, record that fact and abort the loop.
            IF KEY
            THEN
DROP -1 'STAT' STO
63 'i' STO
            END
          NEXT
        \>>
      \>>
@ If the clock display is enabled, draw the time at the bottom of the
@ map (erasing the old time first, since GXOR is being used.)
      IF -40 FS?
      THEN OSTR
XTime DT TM TSTR
DUP 'OSTR' STO
XTime
      END
    \>>


@ Main driver program: Compute the relevant astronomical variables and
@ then draw the map.
  SunClk
@ First make sure the right mathematical modes are in effect.
    \<< RCLF -19 CF
-15 CF -16 CF -17
SF -18 CF CLLCD
@ Display an explanatory screen during the initial computation.
"       SunClock" 1
DISP
"  Figuring widths of"
3 DISP
"  daylight bands."
4 DISP
"[Preliminaries]" 7
DISP PICT PURGE
@ Put the completely dark map into PICT.
PICT { # 2d # 0d }
World REPL ""
@ Note that no time string has yet been displayed.
'OSTR' STO
@ This is the animation loop; if "NOW" was pressed, it will be
@ executed multiple times.
      DO JTime 0 0
0 0 \-> jt gt ra dec
xl
        \<< jt DUP
GMST 'gt' STO 0
@ Figure out the location of the sun at the specified instant.
SunPos DROP DROP
'dec' STO 'ra' STO
@ We only care about the angular information.
          IF dec
@ If the declination has changed more than the threshold value, or if
@ it has changed sign, recompute the table of daylight widths.
OldDec - ABS Thres
\>= dec SIGN OldDec
SIGN \=/ OR
          THEN dec
PrjIll
          END
@ Unless the user has already hit a key (requesting termination), draw
@ the resulting map.
          IF STAT 0
\>=
          THEN {
# 0d # 0d } PVIEW
@ Figure out the longitude of noon.
ra gt 15 * - 180 +
FixAng 126 * 360 /
@ Draw the current map.
FLOOR DUP Draw
@ Keep track of the last location of noon.
'ONOON' STO
          END
        \>>
@ If the user has hit a key, end the animation loop if it was active.
        IF KEY
        THEN DROP
-1 'STAT' STO
        END
@ If we're animating, update the date and time and record the old
@ widths table so they can be erased properly on the next update.
        IF STAT 0 >
        THEN DATE
TIME 'TM' STO 'DT'
STO GTime WTAB
'OTAB' STO
        END
@ If we're animating, loop back and recompute/redraw.
      UNTIL STAT 0
\<=
@ Get rid of extraneous variables.
      END STAT {
OTAB ONOON OSTR
PPAR STAT } PURGE
@ If the user hasn't hit a key but we're ending because "THEN" was the
@ calling program (one-shot date display), freeze the screen so the
@ user can enjoy their map.
      IF 0 \>=
      THEN 7 FREEZE
      END STOF
    \>>


@ Figure out what widths of sunlight fall on the latitude bands of the
@ map. Given the sun's declination in level 1.
  PrjIll
    \<< \-> dec
      \<< 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 \->
i ilon ilat lilon
lilat xt m x y z th
lon lat s c
@ Initialize the width table to "darkness". This was originally an
@ array, but 64 real numbers took too much space, so I changed it to a
@ string. 255 means no light, other numbers are 1/2 the width to
@ illuminate.
        \<< "" 1 64
          FOR c 255
CHR +
          NEXT
@ Flag 6 is set the first time through, so interpolation isn't attempted.
'WTAB' STO 6 SF dec
@ Build the transformation for the declination.
D\->R NEG SIN 's' STO
dec D\->R NEG COS 'c'
STO \pi \->NUM 2 / NEG
'th' STO
@ Increment over a semicircle of illumination.
          DO s NEG
@ Transform the point through the declination rotation.
th SIN * 'x' STO th
COS 'y' STO c th
SIN * 'z' STO y x
@ Transform the resulting coordinate through the map projection to
@ obtain screen coordiantes.
ATAN2 R\->D 'lon' STO
z ASIN R\->D 'lat'
STO 62 lat 90 +
.344444444444 * -
IP 'ilat' STO lon
.35 * IP 'ilon' STO
            IF 6
FS?C
            THEN
@ First time through; just record the previous coordinates.
ilon 'lilon' STO
ilat 'lilat' STO
            ELSE
@ Interpolate between the current and previous points if we've jumped
@ more than one map line down.
              IF
lilat ilat ==
              THEN
WTAB 62 ilat - ilon
1 MAX CHR REPL
'WTAB' STO
              ELSE
ilon lilon - ilat
lilat - / 'm' STO
lilat 'i' STO
WHILE i ilat \=/
REPEAT lilon i
lilat - m * .5 +
FLOOR + 'xt' STO
WTAB 62 i - xt 1
MAX CHR REPL 'WTAB'
STO ilat lilat -
SIGN 'i' STO+
END
              END
ilon 'lilon' STO
ilat 'lilat' STO
            END \pi
@ Back at the outer loop; give the user feedback about how the
@ computation is progressing.
\->NUM Res / 'th'
STO+ "[" th \pi \->NUM
2 / + \pi \->NUM / 100
* 100 MIN FLOOR
\->STR + "%]" + 7
DISP
@ If the user wishes to terminate this computation, record the fact
@ and abort the loop.
            IF KEY
            THEN
DROP 3 'th' STO -1
'STAT' STO 1000
'OldDec' STO
            END
@ Keep computing until the full semicircle is done.
          UNTIL th
\pi \->NUM 2 / .001 + >
          END
@ One edge of the map is fully illuminated, but may have been missed
@ in the above loop, so tweak the widths at that end.
          IF dec 0
<
@ South pole in daylight
          THEN -1
64
          ELSE 1 1
@ North pole in daylight.
          END
'ilat' STO 'lilat'
STO
@ Don't bother if the user's already aborting.
          IF STAT 0
\>=
@ Loop from the edge of the map until a non-dark line is found,
@ illuminating as you go.
          THEN ilat
31
            FOR j
              IF
WTAB j DUP SUB NUM
255 \=/
              THEN
DO WTAB j 63 CHR
REPL 'WTAB' STO
  IF j ilat ==
  THEN lilat 'j'
STO
  END lilat NEG 'j'
STO+
UNTIL j 0 ==
END 31 'j' STO
              END
lilat
            STEP
@ All done; keep track of the declination value used in this
@ computation, so it can be avoided next time if the current
@ declination is close enough.
dec 'OldDec' STO
          END
        \>>
      \>>
    \>>


@ Calculate the position of the sun. Level 2 contains the Julian time
@ of the instant for which the position is desired, and level 1 should
@ be nonzero if the apparent position (corrected for nutation and
@ aberration) is desired.
@ The Sun's longitudinal position (in degrees) is returned in level 4;
@ divide by 15 to get hours. The declination is returned in level 3.
@ The radius vector (in astronomical units) is returned in level 2,
@ and the Sun's longitude (true or apparent, as desired) is returned
@ as degrees in level 1.
  SunPos
    \<< \-> jd ap
      \<< 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 \->
t t2 t3 l m e ea v
th om eps ra dec rv
slong
        \<< jd
2415020 - 36525 /
@ Convert time to Julian centuries of 36525 ephemeris days, measured
@ from the epoch 1900 January 0.5 ET.
DUP 't' STO DUP DUP
* DUP 't2' STO *
't3' STO 279.69668
36000.76892 t * +
.0003025 t2 * +
FixAng 'l' STO
@ 'l' gets the geometric mean longitude of the Sun, referred to the
@ mean equinox of the date. Now compute the sun's mean anomaly.
358.47583
35999.04975 t * +
.00015 t2 * -
.0000033 t3 * -
FixAng 'm' STO
@ Now take the eccentricity of the Earth's orbit into account.
.01675104 .0000418
t * - .000000126 t2
* - 'e' STO m e
@ Compute the eccentric anomaly.
Kepler 'ea' STO 1 e
@ Compute the true anomaly.
+ 1 e - / \v/ ea 2 /
TAN * ATAN R\->D 2 *
FixAng 'v' STO l v
@ Compute the sun's true longitude.
+ m - 'th' STO
@ Obliquity of the ecliptic:
23.452294 .0130125
t * - .00000164 t2
* - .000000503 t3 *
+ 'eps' STO
          IF ap
@ Corrections for sun's apparent longitude, if desired.
          THEN
259.18 1934.142 t *
- FixAng 'om' STO
.00569 .00479 om
D\->R SIN * - NEG
'th' STO+ .00256 om
D\->R COS * 'eps'
STO+
          END th
@ Record sun's longitude and radius vector, then compute coordinates.
'slong' STO
1.0000002 1 e e * -
* 1 e v D\->R COS * +
/ 'rv' STO eps D\->R
COS th D\->R SIN * th
D\->R COS ATAN2 R\->D
FixAng 'ra' STO eps
D\->R SIN th D\->R SIN
* ASIN R\->D 'dec'
STO ra dec rv slong
        \>>
      \>>
    \>>


@ Solve the equation of Kepler.
  Kepler
    \<< \-> m ecc
      \<< 0 0 \-> e d
        \<< m D\->R DUP
'm' STO 'e' STO
          DO e e
SIN ecc * - m - 'd'
STO d 1 e COS ecc *
- / NEG 'e' STO+
          UNTIL d
ABS .000001 \<=
          END e
        \>>
      \>>
    \>>

@ Compute arctan(y/x), with y in level 2 and x in level 1. Returns
@ appropriate quadrant treating y and x as rectangular coordinates
@ being converted to polar coordinates.
  ATAN2
    \<< SWAP \->V2 -16
SF V\-> -16 CF SWAP
DROP
    \>>

@ Translate a big angle back in to the range 0-360 degrees.
  FixAng
    \<< DUP 360 /
FLOOR 360 * -
    \>>

@ Show the currently chosen date and time and prompt the user for
@ their next action (used by the "THEN" submenu).
  ShwTm
    \<< CLLCD DT TM
TSTR 4 DISP
"Choose time, push GO"
3 DISP 2 FREEZE
    \>>
  GMST
    \<< \-> jd
      \<< jd .5 +
FLOOR .5 - 2415020
- 36525 / 0 \-> t th0
        \<< 6.6460656
2400.051262 t * +
.00002581 t t * * +
'th0' STO jd .5 +
DUP FLOOR - 24 *
1.002737908 * th0 +
DUP 24 / FLOOR 24 *
-
        \>>
      \>>
    \>>

@ Convert the chosen date and time to Greenwich time.
  GTime
    \<< RCLF -42 CF
DT 'GD' STO TM TZ -
DUP
      IF 24 \>=
      THEN 24 - GD
1 DATE+ 'GD' STO
      END 'GT' STO
STOF
    \>>

@ Compute the chosen Julian time, as a date and day fraction.
  JTime
    \<< JDate .5 - GT
HMS\-> 24 / +
    \>>

@ Compute the Julian date for the chosen date.
  JDate
    \<< GD DUP IP
SWAP FP 100 * DUP
IP SWAP FP 10000 *
0 \-> m d y c
      \<<
        IF 'm>2'
        THEN -3 'm'
STO+
        ELSE 9 'm'
STO+ -1 'y' STO+
        END y 100 /
IP DUP 'c' STO 100
* NEG 'y' STO+ d c
146097 * 4 / IP + y
1461 * 4 / IP + m
153 * 2 + 5 / IP +
1721119 +
      \>>
    \>>

@ Holds the currently chosen time:
  TM 22.3138814208

@ Holds the currently chosen date:
  DT 11.151991

@ Greenwich time and date for above values:
  GT 3.3138814208
  GD 11.161991

@ Stores the declination for the last computed daylight width table.
  OldDec
-18.5917760031

@ The most recently computed width table, as a string for space
@ reasons (about 70 bytes instead of 500, and there are two of these
@ around while animation is active.)
  WTAB
C$ 64 \255\255\255\255\255\255
   !!!"""##$$$%&'')*+,/5?????????
END